//单纯形算法
//simplex.cpp

//一个线性不等式组是这样的：
//-2*x1 + 8*x2 + 10*x4 >= 50
//5*x1 + 2*x2 >= 100
//3*x1 - 5*x2 + 10*x3 - 2*x4 >= 25
//在这个线性不等式组上给定一个目标：x1 + x2 + x3 + x4最小
//希望在满足线性不等式组的情况下求出令x1+x2+x3+x4最小的一个方案，这类问题称为
//线性规划问题

//1)线性规划的一般形式
//线性规划的问题总是希望最优化一个满足一组线性不等式构成的约束条件的线性函数
//线性函数f定义为：f(x1, x2, ..., xn) = a1x1 + a2x2 + ... + anxn
//其中a1到an是常数，x1到xn是变量
//将等式f(x1, x2, ..., xn) = b称为一个线性等式，其中b是常数
//将这样的不等式f(x1, x2, ..., xn) <= b称为一个线性不等式，其中b是常数
//当未知量x1到xn满足以下线性约束条件：
//f1(x1, x2, ..., xn) <= b1
//f2(x1, x2, ..., xn) <= b2
//......
//fm(x1, x2, ..., xn) <= bm
//然后给出一个目标z = a1x1 + a2x2+ ... + anxn
//其中b1到bm和a1到an都是常数
//希望求出能够使z最大或最小的一组x = (x1, x2, ..., xn)，这一组x的解是一个最优解
//但是在线性约束条件中不能出现单纯的 < 或 > 号
//因为开区间会使得x的解无限趋近于临界值但是永远无法等于，使得最优解失去意义
//所以线性约束条件中只能出现 >=，<= 和 = 号
//在一确定线性约束下，可以求目标z的最小化，也可以求其最大化
//而求两种最优化使用的求解算法是一样的，本节只考虑求第一种的最大化
//
//2)线性规划映射坐标系
//考虑这样一个简单的线性规划问题：
//约束条件：
//4*x1 - x2 <= 8
//2*x1 + x2 <= 10
//5*x1 - 2*x2 >= -2
//x1 >= 0
//x2 >= 0
//目标为z = x1 + x2的最大化
//当我们将x1和x2映射在二维平面坐标系(笛卡尔坐标系)上的x和y坐标时
//我们会发现约束条件实际组成了一个由五条直线划分出的一块平面区域，目标z也是一条直线
//当z取值不同时，z直线可以在坐标系平面内任意移动，但斜率始终不变
//当z直线穿过平面区域，则在平面区域内的线段上的任意一点都满足约束条件
//称为线性约束的一个可行解，当然不一定是所求的目标z
//可以轻松的发现，当x2，即坐标系中的y坐标最大时，目标z具有最大值
//当y坐标(x2)最小时目标z具有最小值
//将z直线向右上方移动，直到z直线与平面区域的最右上角处于邻接位置时，得到z的最大化
//将z直线向左下方移动，直到z直线与平面区域的最左下角处于邻接位置时，得到z的最小化
//当z直线取最大或最小值的临界位置与平面区域的一条直线完全重合
//这时z的最大最小化方案有无穷多，即这两条直线重合的线段上理论上有无穷多个点
//线性规划在二维平面坐标系中映射的平面区域的一个有趣的性质是它一定是一个凸多边形
//而不会出现凹的边，因为组成该平面区域的是直线，而没有射线或线段
//
//3)单纯形(Simplex)
//当线性规划中的未知量个数小于等于2时
//总可以使用平面坐标系的办法直观的看到线性约束映射的平面图形
//这些涉及到解析几何中直线和平面区域的知识，读者可以参考一些高中数学教材
//当未知量的数量是3个时
//线性约束映射为空间中三维坐标系的一块立体区域
//而目标z映射为三维坐标系中的一条直线
//当z直线穿过立体区域时，穿过的线段上所有点都是线性约束的可行解
//而在临界位置处目标z会取得最大值或最小值
//当未知量的数量为n个时
//线性规划问题所映射的空间是一个n维空间坐标系
//具有n个未知量的线性约束在n维空间坐标系中映射的空间区域就是该线性规划的单纯形
//与二维坐标系中的平面区域有相同的性质
//单纯形一定是凸的空间区域，而不会存在凹的部分
//目标z映射的z直线与单纯形的临界位置相交处的顶点就是目标z的最优解
//
//4)单纯形算法
//当单纯形是凸的
//将z直线沿着单纯形的一条边(即一条临界直线)移动到大于等于当前目标z的位置
//这个位置一般是单纯形的一个顶点
//直到z直线所在的顶点是这个局部中的最大值，由于单纯形是凸的，局部最优也是全局最优
//则该顶点也是目标z的最大化，算法结束
//
//5)标准型(Standard)
//标准型是线性规划的一个固定格式：
//线性约束总是满足：
//a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn <= b1
//a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn <= b2
//......
//a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn <= bm
//x1 >= 0
//x2 >= 0
//......
//xn >= 0
//目标为最大化z = c1x1 + c2x2 + ... + cnxn
//用前一节中的矩阵来表示线性不等式时，则标准型可以简单的描述为：
//线性约束：
//Ax <= b
//x >= 0
//目标为最大化z = cx
//其中A是m*n的矩阵，x和c是n维向量，b是m维向量
//一个满足所有线性规划的向量x是可行解，否则是不可行解
//令目标z最大化的解x是最优解，而x取最优解时的目标值z是最优目标值
//若线性规划没有可行解则称该线性规划不可行，否则是可行的
//若线性规划存在无穷多的最优解则称该线性规划是无界的
//
//6)一般形转化为标准型
//一般线性规划中并不一定具备标准型的格式，有四种格式不一样的地方：
//i)目标函数z = c1x1 + c2x2 + ... + cnxn求最小化，而不是最大化
//ii)变量x不具有非负性约束x >= 0
//iii)可能存在等式约束，即存在a1x1 + a2x2 + ... + anxn = b
//iv)可能存在大于不等式约束，即存在a1x1 + a2x2 + ... + anxn >= b
//对于以上四种情况，依次进行如下变换：
//i)将目标函数z的系数取负值，变为z' = -c1x2 - c2x2 - ... - cnxn
//将求z的最小化转化为求目标z'的最大化
//ii)将不满足非负性约束的变量x用量两个新的变量代替，即x = x' - x''
//并且增加非负性条件x' >= 0和x'' >= 0，将原线性规划中的x都由x'-x''代替
//iii)将等式约束转化为等效的两个小于等于的不等式约束
//将原来的a1x1 + a2x2 + ... + anxn = b改为以下两个小于等于不等式约束
//a1x1 + a2x2 + ... + anxn <= b和
//-a1x1 - a2x2 - ... - anxn <= -b，即a1x1 + a2x2 + ... + anxn >= b的小于等于形式
//iv)将大于等于不等式取负，转化为小于等于的形式
//将原来的a1x1 + a2x2 + ... + anxn >= b转化为
//-a1x1 - a2x2 - ... - anxn <= -b
//通过上面四种转化可以将所有线性不等式都转化为标准型
//
//7)松弛形(Relaxed)
//松弛形是一组全部由等式约束描述的线性规划，它可以由标准型转化而来
//对于这样的一个标准型：
//a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn <= b1
//a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn <= b2
//......
//a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn <= bm
//x1 >= 0
//x2 >= 0
//......
//xn >= 0
//目标为最大化z = c1x1 + c2x2 + ... + cnxn
//为标准型添加m个新的未知量，依次为x1',x2',...,xm'，使标准型转化为松弛形：
//x1' = b1 - a(1,1)x1 - a(1,2)x2 - ... - a(1,n)xn
//x2' = b2 - a(2,1)x1 - a(2,2)x2 - ... - a(2,n)xn
//...
//xm' = bm - a(m,1)x1 - a(m,2)x2 - ... - a(m,n)xn
//x1 >= 0
//......
//xn >= 0
//x1' >= 0
//......
//xm' >= 0
//(目标)z = 0 + c1x1 + c2x2 + ... + cnxn，其中0对应等式中的系数b
//所有未知量都受到非负性约束，使得松弛形与标准型是等价的线性不等式组
//等式左边新添加的未知量x'称为基本变量，基本变量的集合为B(Basic)
//等式右边的未知量x称为非基本变量，非基本变量的集合为N(Non-Basic)
//也可以写作矩阵形式：
//x' = b - Ax
//x >= 0, x' >= 0
//其中x'是m维向量
//
//8)线性规划的应用
//线性规划可以用来解决很多问题，包括图论一章中的最短路径，最大流等问题
//最短路径可以看作已知图中各节点的距离关系，求最小化目标距离函数的问题
//最大流可以看作已知流网络中各节点的容量关系，求最大化目标最大流的问题
//
//9)单纯形算法
//让我们回头再来看单纯形算法
//单纯形是有n个未知量的线性规划的标准型在n维空间中映射的空间区域
//它总是一个凸的空间区域，而不存在凹的部分
//单纯形算法的基本思想是将线性目标函数z沿着单纯形的边移动，找出局部最大化的目标函数z
//此局部最大的目标z即为单纯形全局最大的目标z
//而单纯形算法是通过在松弛形上进行一个循环的变换来实现的
//这里本文讲解一个引用自算法导论中单纯形算法的例子
//考虑这样一个格式为松弛形的线性规划：
//最大化z = 3*x1 + x2 + 2*x3						(1)
//x4 = 30 - x1 - x2 - 3*x3							(2)
//x5 = 24 - 2*x1 - 2*x2 - 5*x3						(3)
//x6 = 36 - 4*x1 - x2 - 2*x3						(4)
//以上的所有未知量x>=0，等号左边的是基本解，等号右边的是非基本解
//现在开始重复的对该松弛形重复以下步骤直到目标z无法继续增大
//令所有右边的非基本解(x1, x2, x3)为0
//此时z = 0，(x1, x2, x3, x4, x5, x6) = (0, 0, 0, 30, 24, 36)
//第一次操作：
//在当前式(1)中选择能使z最快增大的未知量，即系数为正且值最大的未知量
//这个未知量称为主元，选取的过程称为选主元
//选择x1，在式(2)(3)(4)中计算x1最大能取到的值
//由于x4，x5，x6非负，所以x1 = 9为最大值，在式(4)处取得，即式(4)处x1最紧
//故将x6和x1交换，将x1放在等号左边作为基本解，x6放在等号右边作为非基本解
//式(4)经过等价变换得：x1 = 9 - x2/4 - x3/2 - x6/4，再代入式(1)(2)(3)，得到新松弛形：
//z = 27 + x2/4 + x3/2 - 3*x6/4						(5)
//x1 = 9 - x2/4 - x3/2 - x6/4						(6)
//x4 = 21 - 3*x2/4 - 5*x3/2 + x6/4					(7)
//x5 = 6 - 3*x2/2 - 4*x3 + x6/2						(8)
//令所有右边的非基本解(x2, x3, x6)为0
//此时z = 27(z增大了)，(x1, x2, x3, x4, x5, x6) = (9, 0, 0, 21, 6, 0)
//第二次操作：
//在当前式(5)中选择能使z最快增大的未知量，注意x6的系数是负值且x6非负所以不能选x6
//选择x3，在式(6)(7)(8)中计算x3最大能取到的值
//由于x1，x4，x5非负，所以x3 = 1.5为最大值，在式(8)处取得，即式(8)处x3最紧
//将x3与x5交换位置，x3放在左边x5放在右边，再带入式(5)(6)(7)(8)中可得新松弛形：
//z = 111/4 + x2/16 - x5/8 - 11*x6/16				(9)
//x1 = 33/4 - x2/16 + x5/8 - 5*x6/16				(10)
//x3 = 3/2 - 3*x2/8 - x5/4 + x6/8					(11)
//x4 = 69/4 + 3*x2/16 + 5*x5/8 - x6/16				(12)
//令所有右边的非基本解(x2, x5, x6)为0
//此时z = 111/4，(x1, x2, x3, x4, x5, x6) = (33/4, 0, 3/2, 69/4, 0, 0)
//第三次操作：
//在当前式(9)中选择能使z最快增大的未知量
//事实上这时只能选择x2了，因为其他两未知量的系数为负值
//x2在式(10)处最大可取132，在式(11)处最大可取4
//而在式(12)处最大可取无穷大，因为式(12)中x2系数为正，x4随着x2增大而增大
//这时x2选择4，即x2在式(11)处最紧
//交换x2与x4的位置，并且变换式(9)(10)(11)(12)，得到新松弛形：
//z = 28 - x3/6 - x5/6 - 2*x6/3						(13)
//x1 = 8 + x3/6 + x5/6 - x6/3						(14)
//x2 = 4 - 8*x3/3 - x*x5/3 + x6/3					(15)
//x3 = 18 - x3/2 + x5/2								(16)
//令所有右边的非基本解(x3, x5, x6)为0
//此时z = 28，(x1, x2, x3, x4, x5, x6) = (8, 4, 0, 18, 0, 0)
//式(13)中所有未知量的系数都是负数，且所有未知量都是非负的
//也就是说，只有当非基本解都取0时，z有最大值，即目标z的最大值即为28，算法结束
//带回原始的目标函数z = 3*x1 + x2 + 2*x3有x1 = 8，x2 = 4，x3 = 0时z = 28
//
//当然线性规划有解的情况中也分有唯一解和有无穷多解，单纯形算法无法找出无穷多解
//而只能给出无穷多解中的一个
//
//10)线性规划的无界
//在单纯形算法重复进行的时候，如果有一次的步骤中出现这样的情况
//在z中选择出一个未知量x，在当前松弛形的等式中查找那个能让x最紧的等式
//若无法找出一个能让x最紧的等式，即这个x在所有等式中都能取得无穷大
//则该线性规划的是无界的，即目标函数z向上可以取得无穷大，算法结束
//
//11)线性规划的无解
//判定线性规划的无解需要构造一个辅助线性规划，可以将它轻松的转化为松弛形
//考虑一个标准型L：
//线性约束：
//a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn <= b1
//a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn <= b2
//......
//a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn <= bm
//x1 >= 0
//x2 >= 0
//......
//xn >= 0
//目标为最大化z = c1x1 + c2x2 + ... + cnxn
//
//向标准型L中添加未知量x0构造为对应的辅助线性规划Laux：
//线性约束：
//a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn - x0 <= b1
//a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn - x0 <= b2
//......
//a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn - x0 <= bm
//x0 >= 0
//x1 >= 0
//x2 >= 0
//......
//xn >= 0
//目标为最大化z = -x0
//将辅助线性规划构造为松弛形Raux，再预先进行一次交换操作，将未知量x0换到等号左边
//交换的行是系数b最小的那行，然后对这个松弛形Raux用单纯形算法求解
//当且仅当Laux的最有目标值z为0时，L是可行的，否则L无解
//
//12)初始化
//一般格式的线性规划可以简单的转化为标准型，因此本文将标准型作为线性规划的输入
//标准型无法直接应用单纯形算法，而需要转化为松弛形
//在转化之前先构造辅助线性规划进行一次单纯形算法，判断出该线性规划是否无解
//若有解则应用单纯形算法求出最优方案
//
//我忠实的实现了算法导论第29章线性规划中的关于单纯形算法的描述
//注意代码中各数组的下标都是从0开始的，所有的下标都比上面讲解中的小1
//本节最后的测试用例将无解，无界，有唯一解和有无穷多解的情况都进行了测试

#include "general_head.h"
#include "matrix.h"
bool initialize_l(matrix& e, int& m, int& n, double *c);
void construct_aux(matrix& e, int m, int& n, double *c);
void standard2relaxed(matrix& e, int m, int n, double *c);
bool simplex_core(matrix e, int m, int n, double *c, double *x, double& z);
bool select_pivot(double *c, int *non_basic, int mn, int& pivot);
void rewrite_matrix(matrix& e, int m, int n, int pivot, int row,
		double *c, int *basic, int *non_basic);

int simplex(matrix e, int m, int n, double *c, double *x, double& z)
{//矩阵e是m行n+1列的，0列到n-1列是不等号左边的矩阵A，最后一列n列是m维向量的系数b
 //c是目标函数z的n维向量
 //该线性规划是标准型，所有未知量x都是非负的
 //对该线性规划求解，无解返回-1，无界返回1
 //若有解返回0，返回最优解z，最优方案存储于数组x中
	if(!initialize_l(e, m, n, c))
		return(-1);
	if(!simplex_core(e, m, n, c, x, z))
		return(1);
	else
		return(0);
}
bool initialize_l(matrix& e, int& m, int& n, double *c)
{//矩阵e是m行n+1列的，0列到n-1列是不等号左边的矩阵A，最后一列n列是m维向量的系数b
 //c是目标函数z的n维向量
 //该线性规划是标准型，所有未知量x都是非负的
 //返回标准型是否有解，若有解则将矩阵A构造为松弛形
	//构造出标准型L的松弛形R，对松弛形R进行预测试
	//若松弛形R的系数向量b都大于等于0，则该线性规划至少存在x全为0的解
	//直接将标准型L转化为松弛形R，返回
	matrix r(e);
	int r_m(m), r_n(n);
	double r_c[MAX];
	standard2relaxed(r, r_m, r_n, r_c);
	bool t_flag(true);
	for(int i = 0; i < r_m; ++ i)
		if(r.m_m[i][r_m + r_n] < 0)
			t_flag = false;
	if(t_flag){
		standard2relaxed(e, m, n, c);
		return(true);
	}

	matrix aux(e);
	int m_aux(m), n_aux(n);
	double c_aux[MAX], x_aux[MAX], z_aux;
	//将标准型L构造为对应的辅助线性代数的标准型Laux
	construct_aux(aux, m_aux, n_aux, c_aux);
	//将标准型Laux构造为松弛形
	standard2relaxed(aux, m_aux, n_aux, c_aux);

	//将松弛形的未知量x0与系数b最小的那行的基本解交换位置
	double tight(INF), row(0);
	for(int i = 0; i < m; ++ i)
		if(tight > aux.m_m[i][m_aux + n_aux]){
			tight = aux.m_m[i][m_aux + n_aux];
			row = i;
		}
	int basic[MAX], non_basic[MAX];
	memset(basic, -1, MAX * sizeof(int));
	memset(non_basic, -1, MAX * sizeof(int));
	//下标n到m+n-1(共m个)是基本解，i行等号左边基本解的下标是i+n
	for(int i = 0; i < m_aux; ++ i)
		basic[i] = i + n_aux;
	//初始时下标0到n-1是非基本解
	for(int i = 0; i < n_aux; ++ i)
		non_basic[i] = 1;
	rewrite_matrix(aux, m_aux, n_aux, n_aux - 1, row, c_aux, basic, non_basic);

	//对Laux松弛形应用单纯形算法
	simplex_core(aux, m_aux, n_aux, c_aux, x_aux, z_aux);
	//若Laux松弛形的最优解z_aux != 0，返回无解
	if(z_aux != 0)
		return(false);
	//若Laux松弛形的最优解z_aux = 0，返回有解以及标准型L的松弛形
	else{
		standard2relaxed(e, m, n, c);
		return(true);
	}
}
void construct_aux(matrix& e, int m, int& n, double *c)
{//矩阵e是m行n+1列的，0列到n-1列是不等号左边的矩阵A，最后一列n列是m维向量的系数b
 //c是目标函数z的n维向量
 //将标准型L构造为对应的辅助线性代数的标准型Laux
 //标准型L：
 //a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn <= b1
 //a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn <= b2
 //......
 //a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn <= bm
 //目标为最大化z = c1x1 + c2x2 + ... + cnxn，所有未知量都非负
 //添加未知量xn+1，而不是x0，可以直接添加在数据结构的最后，而不需要移动所有元素
 //得到对应的辅助线性代数Laux：
 //a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn - xn+1 <= b1
 //a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn - xn+1 <= b2
 //......
 //a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn - xn+1 <= bm
 //目标为最大化z = -xn+1，所有未知量都非负
	//注意代码中各数组的下标都是从0开始的，所有的下标都比讲解中的减1
	//设置目标函数z的系数向量c
	memset(c, 0, MAX * sizeof(int));
	c[n] = -1;
	++ n;
	for(int i = 0; i < m; ++ i){
		//将每一行的系数b向后移动一个位置
		e.m_m[i][n] = e.m_m[i][n - 1];
		//系数b的原位置变为新未知量xn+1的位置
		e.m_m[i][n - 1] = -1;
	}
}
void standard2relaxed(matrix& e, int m, int n, double *c)
{//矩阵e是m行n+1列的，0列到n-1列是不等号左边的矩阵A，最后一列n列是m维向量的系数b
 //c是目标函数z的n维向量
 //标准型:
 //a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn <= b1
 //a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn <= b2
 //......
 //a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn <= bm
 //目标为最大化z = c1x1 + c2x2 + ... + cnxn，所有未知量都非负
 //添加m个未知量，将标准型构造对应的松弛形：
 //z = c1x1 + c2x2 + ... + cnxn
 //xn+1 = b1 - a(1,1)x1 - a(1,2)x2 - ... - a(1,n)xn
 //xn+2 = b2 - a(2,1)x1 - a(2,2)x2 - ... - a(2,n)xn
 //......
 //xn+m = bm - a(m,1)x1 - a(m,2)x2 - ... - a(m,n)xn
 //所有未知量都非负
	//为了方便我并未改变矩阵A中的未知量的系数的正负号
	//只需将系数向量b移动到n+m列，前面0列到m+n-1列是为m+n个未知量腾出位置
	for(int i = 0; i < m; ++ i){
		e.m_m[i][m + n] = e.m_m[i][n];
		e.m_m[i][n] = 0;
	}
	e.m_col = m + n;
}
bool simplex_core(matrix e, int m, int n, double *c, double *x, double& z)
{//矩阵e是m行n+1列的，0列到n-1列是不等号左边的矩阵A，最后一列n列是m维向量的系数b
 //c是目标函数z的n维向量
 //输入的线性规划是松弛形：
 //z = c1x1 + c2x2 + ... + cnxn
 //xn+1 = b1 - a(1,1)x1 - a(1,2)x2 - ... - a(1,n)xn
 //xn+2 = b2 - a(2,1)x1 - a(2,2)x2 - ... - a(2,n)xn
 //......
 //xn+m = bm - a(m,1)x1 - a(m,2)x2 - ... - a(m,n)xn
 //等号左边是基本解，等号右边是非基本解，所有未知量都非负
 //判断该松弛形是否有最优解，若有解返回最优解z，将最优方案存储于数组x中
	//注意代码中各数组的下标都是从0开始的，所有的下标都比讲解中的减1
	//basic[i] = j指代矩阵中i行等式左边基本解的下标为j，数组从0开始故用-1表示没有使用
	//non_basic[i]指代未知量xi+1是非基本解
	int basic[MAX], non_basic[MAX];
	memset(basic, -1, MAX * sizeof(int));
	memset(non_basic, -1, MAX * sizeof(int));
	//下标n到m+n-1(共m个)是基本解，i行等号左边基本解的下标是i+n
	for(int i = 0; i < m; ++ i)
		basic[i] = i + n;
	//初始时下标0到n-1是非基本解
	for(int i = 0; i < n; ++ i)
		non_basic[i] = 1;
	//初始化解x
	memset(x, 0, MAX * sizeof(double));
	for(int i = 0; i < m; ++ i)
		x[basic[i]] = e.m_m[i][m + n];
	int pivot;
	while(select_pivot(c, non_basic, m + n, pivot)){
		//若当前目标函数z中非基本解的系数还有正数，则返回最大正系数的主元下标pivot
		//若所有未知量系数都是负数，则无法找到更大的z，算法结束
		double tight(INF), row(0);
		for(int i = 0; i < m; ++ i){
			//求每一行中pivot主元的最大取值，找出所有最大值中最紧的一行row
			if(e.m_m[i][pivot] <= 0)
				//跳过非正数的pivot，此处主元与等号左边的基本解正负号相同
				//即此处主元可以取无穷大INF
				continue;
			double tmp = e.m_m[i][m + n] / e.m_m[i][pivot];
			if(tight > tmp){
				tight = tmp;
				row = i;
			}
		}
		if(tight == INF)
			return(false);
		//对当前矩阵根据主元pivot进行重写
		rewrite_matrix(e, m, n, pivot, row, c, basic, non_basic);
		//重写解x
		memset(x, 0, MAX * sizeof(double));
		for(int i = 0; i < m; ++ i)
			x[basic[i]] = e.m_m[i][m + n];
	}
	z = c[m + n];
	return(true);
}
bool select_pivot(double *c, int *non_basic, int mn, int& pivot)
{//选择目标函数z中系数最大的未知量作为主元
 //若所有系数都小于等于0则没有主元
	double tmp(0.0);
	bool p_flag(false);
	for(int i = 0; i < mn; ++ i)
		if(non_basic[i] == 1 && tmp < c[i]){
			tmp = c[i];
			pivot = i;
			p_flag = true;
		}
	if(p_flag)
		return(true);
	else
		return(false);
}
void rewrite_matrix(matrix& e, int m, int n, int pivot, int row,
		double *c, int *basic, int *non_basic)
{//对松弛形进行重写：
 //z = c1x1 + c2x2 + ... + cnxn
 //xn+1 = b1 - a(1,1)x1 - a(1,2)x2 - ... - a(1,n)xn
 //xn+2 = b2 - a(2,1)x1 - a(2,2)x2 - ... - a(2,n)xn
 //......
 //xn+m = bm - a(m,1)x1 - a(m,2)x2 - ... - a(m,n)xn
 //等号左边是基本解，等号右边是非基本解，所有未知量都非负
	//改写主元交换位置的row行
	for(int i = 0; i <= m + n; ++ i)
		if(i != pivot)
			e.m_m[row][i] /= e.m_m[row][pivot];
	e.m_m[row][basic[row]] = 1 / e.m_m[row][pivot];
	e.m_m[row][pivot] = 0;
	//对basic和non_basic数组上的信息标记
	non_basic[pivot] = -1;
	non_basic[basic[row]] = 1;
	basic[row] = pivot;
	//通过row行改写其他每一行
	for(int i = 0; i < m; ++ i)
		if(i != row){
			for(int j = 0; j < m + n; ++ j)
				//将在i行中的主元pivot由改写后的row行中的非基本解来表示
				if(non_basic[j] == 1)
					e.m_m[i][j] -= e.m_m[i][pivot] * e.m_m[row][j];
			e.m_m[i][m + n] -= e.m_m[i][pivot] * e.m_m[row][m + n];
			e.m_m[i][pivot] = 0;
		}
	//将目标函数z也改写
	for(int i = 0; i < m + n; ++ i)
		if(non_basic[i] == 1)
			c[i] -= c[pivot] * e.m_m[row][i];
	c[m + n] += c[pivot] * e.m_m[row][m + n];
	c[pivot] = 0;
}
